A basic introduction to runoutSIM

Introduction

The runoutSIM package is developed for regional-scale runout simulations (spatial distribution and velocity) using random walks. This vignette provides an introduction to applying runoutSIM for regional debris flow simulation using a subcatchment of the upper Maipo river basin, Chile as an example. In this vignette, we will:

Load and visualize spatial data
Simulate runout from single and multiple source cells
Visualize simulation results as rasters and interactive maps

Loading packages and reading spatial data

Spatial handling in the runoutSIM package is supported by the sf (vector data) and terra (gridded/raster data) packages. In this example, we load/read a 12.5 m spatial resolution DEM, debris flow runout polygons, and corresponding source points mapped for each polygon. We also create a hillshade model to help with visualization of the results.

Tip: Currently, this package is designed to work with data in a UTM coordinate reference system (CRS). You will need to project your own data to a local UTM CRS for it to work.

# load packages
library(runoutSIM)
library(terra)
#> Warning: package 'terra' was built under R version 4.4.3
#> terra 1.8.54
library(sf)
#> Warning: package 'sf' was built under R version 4.4.3
#> Linking to GEOS 3.13.0, GDAL 3.10.1, PROJ 9.5.1; sf_use_s2() is TRUE

# Load digital elevation model (DEM)
dem <- rast("C:/GitProjects/runoutSIM/Data/elev_fillsinks_WangLiu.tif")

# Compute hillshade for visualization 
slope <- terrain(dem, "slope", unit="radians")
aspect <- terrain(dem, "aspect", unit="radians")
hill <- shade(slope, aspect, 40, 270)

# Load debris flow runout source points and polygons
source_points <- st_read("C:/GitProjects/runoutSIM/Data/debris_flow_source_points.shp")
#> Reading layer `debris_flow_source_points' from data source 
#>   `C:\GitProjects\runoutSIM\Data\debris_flow_source_points.shp' 
#>   using driver `ESRI Shapefile'
#> Simple feature collection with 73 features and 1 field
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: 389175.6 ymin: 6293926 xmax: 398788.6 ymax: 6327439
#> Projected CRS: WGS 84 / UTM zone 19S
runout_polygons <- st_read("C:/GitProjects/runoutSIM/Data/debris_flow_runout_polygons.shp")
#> Reading layer `debris_flow_runout_polygons' from data source 
#>   `C:\GitProjects\runoutSIM\Data\debris_flow_runout_polygons.shp' 
#>   using driver `ESRI Shapefile'
#> Simple feature collection with 73 features and 1 field
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 389139.1 ymin: 6293864 xmax: 398852.9 ymax: 6327455
#> Projected CRS: WGS 84 / UTM zone 19S

# Load basin boundary
bnd_catchment <- st_read("C:/GitProjects/runoutSIM/Data/basin_rio_olivares.shp")
#> Reading layer `basin_rio_olivares' from data source 
#>   `C:\GitProjects\runoutSIM\Data\basin_rio_olivares.shp' using driver `ESRI Shapefile'
#> Simple feature collection with 1 feature and 15 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 386042.8 ymin: 6293348 xmax: 404405.3 ymax: 6331286
#> Projected CRS: WGS 84 / UTM zone 19S

Creating interactive maps

Using the leafmap() function, we can quickly create an interactive leaflet map for viewing the data. We use the %>% operator to pipe together multiple leafmap() calls to create a combined leaflet map with multiple layers.

Tip: leafmap() has the option to interactively query raster values with a mouse hover. This option is turned on by default, but can result in a considerable increase in the file size of the leaflet widget (e.g. when exported for sharing). It is recommended to have add_image_query = FALSE when saving / exporting the leaflet map.

map <- leafmap(bnd_catchment, color = '#f7f9f9', fill_color = '#FF000000', weight = 4) %>%
    leafmap(runout_polygons) %>%
    leafmap(source_points, color = '#e74c3c') %>%
    leafmap(hill, palette = grey(0:100/100), opacity = 1, add_legend = FALSE, 
            add_image_query = FALSE) %>%
    leafmap(dem, palette = viridis::mako(10), opacity = 0.6, add_image_query = FALSE)
map

Simulate runout from a source cell

Runout paths are simulated from individual source cells, defined as x, y coordinates. For example, if we are simulating runout from a single source cell defined by a source point (vector data), we use sf::st_coordinates() to extract the x, y values.

# Select a single debris flow and source point for the example
runout_polygon <- runout_polygons[31,]

# Get corresponding source point
source_point  <- st_filter(st_as_sf(source_points), st_as_sf(runout_polygon))

# Get coordinates of source point
source_crds <- st_coordinates(source_point)
print(source_crds)
#>             X       Y
#> [1,] 395698.1 6307421

Now that we have the coordinates, we can run the random walk runout simulations using runoutSim(). Where:

mu and md are the sliding friction coefficient and mass-to-drag parameters of the two-parameter friction model controlling runout distance
slp_thresh, exp_div, and per_fct are parameters controlling the dispersion of the random walks
walks is the number of random walk iterations

sim_paths = runoutSim(dem = dem, 
                      xy = source_crds, 
                      mu = 0.08, 
                      md = 140, 
                      slp_thresh = 35, 
                      exp_div = 2.5, 
                      per_fct = 1.95, 
                      walks = 1000)

# plot structure of sim_paths
str(sim_paths)
#> List of 4
#>  $ start_cell    : num 2803625
#>  $ cell_trav_freq: 'table' int [1:1994(1d)] 3 2 3 11 8 5 3 30 31 24 ...
#>   ..- attr(*, "dimnames")=List of 1
#>   .. ..$ all_indices: chr [1:1994] "2774166" "2774167" "2775633" "2775634" ...
#>  $ cell_max_vel  : num [1:1994] 15.2 16.5 11.4 14.3 16 ...
#>  $ prob_connect  : NULL

runoutSim() returns a list containing:

The cell number of the DEM indicating the location of the source cell
Traverse frequencies for each grid cell
Maximum velocities (m/s) for each grid cell

These values are stored as cell numbers (e.g., $all_indices). Although not calculated here, prob_connect provides the probability the source cell connected to a feature of interest.

Converting results to a raster

To visualize the results, we convert the output from runoutSim() to a raster.

## Convert paths to raster with cell transition frequencies
paths_raster <- walksToRaster(sim_paths, dem, method = "freq")

# Plot results
paths_raster <- crop(paths_raster, ext(runout_polygon)+500)
plot(crop(hill, ext(runout_polygon)+500), col=grey(0:100/100), legend=FALSE,
     mar=c(2,2,1,4), main = "Runout traverse frequency")
plot(paths_raster, legend = T, alpha = 0.5, add = TRUE)
plot(st_geometry(runout_polygon), add = TRUE)
plot(source_point, add = TRUE)

In walksToRaster(), we use method = "freq" to return traverse frequencies. We can also return traverse (ECDF) probabilities by using method = "cdf_prob".

We use velocityToRaster() to get a raster of velocity values (m/s).

## Convert paths to raster with max. velocities
vel_raster <- velocityToRaster(sim_paths, dem)

# Plot results
vel_raster <- crop(vel_raster, ext(runout_polygon)+500)
plot(crop(hill, ext(runout_polygon)+500), col=grey(0:100/100), legend=FALSE,
     mar=c(2,2,1,4), main = "Runout velocity")
plot(vel_raster, col = map.pal('plasma'), legend = TRUE, alpha = 0.5, add = TRUE)
plot(st_geometry(runout_polygon), add = TRUE)
plot(source_point, add = TRUE)

Computing connectivity probability

The runoutSim() function also allows us to calculate the connectivity probability of a source cell to feature of interest - i.e. the proportion of walks that intersect a feature from a given source cell.

In this example, we will determine the connectivity probability of a simulated debris flow to a river channel. Using makeConnFeature() function, we will create a matrix indicating with a value of 1 which cells contain the river channel.

# Load river channel polygon (vector)
river_channel <- st_read("C:/GitProjects/runoutSIM/Data/river_channel.shp")
#> Reading layer `river_channel' from data source 
#>   `C:\GitProjects\runoutSIM\Data\river_channel.shp' using driver `ESRI Shapefile'
#> Simple feature collection with 1 feature and 1 field
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 393504.8 ymin: 6293330 xmax: 396806.5 ymax: 6328149
#> Projected CRS: WGS 84 / UTM zone 19S

# Create a connectivity feature for runoutSim that matches the input DEM
feature_mask <- makeConnFeature(river_channel, dem)

Now that we have a connectivity feature, we can re-run runoutSim() with some additional parameter options: source_connect = TRUE allowing to calculate a connectivity probability and connect_feature = feature_mask.

sim_paths = runoutSim(dem = dem, 
                      xy = source_crds, 
                      mu = 0.08, 
                      md = 140, 
                      slp_thresh = 35, 
                      exp_div = 2.5, 
                      per_fct = 1.95, 
                      walks = 1000,
                      source_connect = TRUE,
                      connect_feature = feature_mask)

# Plot structure of sim_paths
str(sim_paths)
#> List of 4
#>  $ start_cell    : num 2803625
#>  $ cell_trav_freq: 'table' int [1:2009(1d)] 3 3 20 14 6 1 41 41 33 20 ...
#>   ..- attr(*, "dimnames")=List of 1
#>   .. ..$ all_indices: chr [1:2009] "2774166" "2775633" "2775634" "2775635" ...
#>  $ cell_max_vel  : num [1:2009] 15 11.7 13.9 16.3 16.7 ...
#>  $ prob_connect  : num 0.883

# Get connectivity probability
sim_paths$prob_connect
#> [1] 0.883

tp_raster <- walksToRaster(sim_paths, dem, method = "cdf_prob")
conn_raster <- connToRaster(sim_paths, dem)

# Crop results
tp_crop <- crop(tp_raster, ext(runout_polygon)+ 500)
conn_crop <- crop(conn_raster, ext(runout_polygon)+ 500)

# Plot results
par(mfrow = c(2,1))
plot(crop(hill, ext(runout_polygon)+500), col=grey(0:100/100), legend=FALSE,
     mar=c(2,2,1,4), main = "Traverse prob.")
plot(st_geometry(river_channel), col = "lightblue", add = TRUE)
plot(tp_crop, legend = TRUE, alpha = 0.5, add = TRUE)
plot(st_geometry(runout_polygon), add = TRUE)
plot(source_point, add = TRUE)


plot(crop(hill, ext(runout_polygon)+500), col=grey(0:100/100), legend=FALSE,
     mar=c(2,2,1,4), main = "Connectivity prob.")
plot(st_geometry(river_channel), col = "lightblue", add = TRUE)
plot(conn_crop, legend = TRUE, alpha = 0.5, add = TRUE)
plot(st_geometry(runout_polygon), add = TRUE)
plot(source_point, add = TRUE)

Simulate runout from multiple-sources

To simulate runout from multiple source cell locations, we run runoutSim() iteratively, changing the xy coordinates during each iteration.

To create a list of xy coordinates from source points (vector), we use the makeSourceList() function.

# Get coordinates of source points to create a source list object
source_l <- makeSourceList(source_points)

# Perform random walk simulation for each source point
sim_runs <- list()

for(i in 1:length(source_l)){
  
  sim_runs[[i]] <- runoutSim(dem = dem, xy = source_l[[i]], 
                              mu = 0.08, 
                              md = 40, 
                              slp_thresh = 40, 
                              exp_div = 3, 
                              per_fct = 1.9, 
                              walks = 1000)
}

We can now convert this list of multiple random walks from different source cells to raster data using walksToRaster() and velocityToRaster().

# Coerce results to a raster
trav_freq <- walksToRaster(sim_runs, method = "freq", dem)
vel_ms <- velocityToRaster(sim_runs, dem, method = "max")
trav_prob <- walksToRaster(sim_runs, method = "max_cdf_prob", dem)


# Plot random walks from mulitple source cells
par(mfrow = c(1,3))
plot(hill, col=grey(0:100/100), legend=FALSE, mar=c(2,2,1,4), 
     main = "Traverse frequency")
plot(trav_freq, add = TRUE, alpha = 0.7)

plot(hill, col=grey(0:100/100), legend=FALSE, mar=c(2,2,1,4), 
     main = "Traverse probability (max.)")
plot(trav_prob, add = TRUE, alpha = 0.7)

plot(hill, col=grey(0:100/100), legend=FALSE, mar=c(2,2,1,4), 
     main = "Velocity (m/s)")
plot(vel_ms, col = map.pal("plasma"), add = TRUE, alpha = 0.7)

We can also view the results through an interactive leaflet map.


# Create interactive leaflet map
sim_map <- 
  # Debris flow observations
  leafmap(runout_polygons, opacity = 0.3, label = "Runout polygons") %>%
  leafmap(source_points, color = '#e74c3c', label = "Source points") %>%
  # Terrain data
  leafmap(hill, palette = grey(0:100/100), opacity = 1, add_legend = FALSE, 
            add_image_query = FALSE, label = "Hillshade") %>%
  leafmap(dem, palette = viridis::mako(10), opacity = 0.6, 
          add_image_query = FALSE, label = "DEM") %>%
  # Modelling results
  leafmap(trav_prob,  palette = viridis::viridis(10, direction = -1), 
          label = "Traverse prob", add_image_query = FALSE) %>%
  leafmap(vel_ms,  palette = viridis::plasma(10, direction = -1),
          label = "Velocity", add_image_query = FALSE)

# Start with mapped zoomed in
sim_map <- leaflet::setView(sim_map, lng = -70.13694, lat = -33.3703, zoom = 13)


# Hide layers
sim_map <- leaflet::hideGroup(sim_map, c("Velocity", "Source points", "DEM", "Hillshade"))

sim_map

Jason Goetz